package scu.maqiang.possion;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.*;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

/**
 * 周期边界下二维Possion方程求解
 */
public class PoissonPeriodic2DT3 {

	public static void main(String[] args) {
		Mesh2T3 mesh2D = new Mesh2T3().square2D(100, 100, T3Type.Cross2);
		mesh2D.scale(Math.PI * 2, Math.PI * 2);
		ScalarFunc sourceFunc = (xy, label, param) -> Math.sin(xy[0] + 0.25 * Math.PI) * Math.cos(xy[1] + 0.25 * Math.PI);
		FES2T31 fs = new FES2T31(mesh2D);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		double[] RHS = new double[fs.GetNdof()];
		fs.assembleStiff(new double[] {1.0}, BVPType.COMMON, A);
		fs.assembleSource(new ScalarFunc[] {sourceFunc}, null, BVPType.COMMON, RHS);
		double[] x = new double[fs.GetNdof()];
		//A.sort();
		//System.out.println(A.display());
		fs.applyBC_Periodic(A, RHS, Direct.Y, 1, 3);
		fs.applyBC_Periodic(A, RHS, Direct.X, 2, 4);
		//fs.applyBC_Periodic(A, RHS, new int[][] {{2, 4}, {1, 3}}, new Direct[] {Direct.X, Direct.Y});
//		System.out.println(A.display());
		//fs.applyBC_MBN(A, RHS, Direct.All, 0.0, 1, 3);
		//A.sort();
		//System.out.println(A.display());
		//System.out.println(A.isSymmetric());
		NewIterSSolver solver = new NewIterSSolver(A);
		solver.PCGSSOR(RHS, x, 1.5, 1);
//		DirectSSolver solver = new DirectSSolver(A, fs.getnBand());
//		NewIterSSolver solver = new NewIterSSolver(A);
//		solver.PLUSolve(RHS, x, 1);
		mesh2D.toTecplot("PeriodicSolution.dat", x);

//		PeriodicCase2();
	}
	
	public static void PeriodicCase2() {
		Mesh2T3 mesh2D = new Mesh2T3().readMesh("mesh2D\\lingXing.msh", MeshFormat.msh);
//		int[][] idxPair = mesh2D.extractPeriodicBoundary(1, 3, new double[] {1.0, -1.0});
		//int[][] idxPair2 = mesh2D.extractPeriodicBoundary(2, 4, new double[] {1.0, 1.0});
		//mesh2D.displayNodes();
		//mesh2D.displayBoundarys();
//		System.out.println(MVO.toString(idxPair));
//		System.out.println(MVO.toString(idxPair2));
				
		FES2T31 fs = new FES2T31(mesh2D);
		
		ScalarFunc sourceFunc = (xy, label, param) -> (xy[0] + xy[1] + 1)*(xy[0] + xy[1] - 1) * (xy[1] - xy[0] + 1) *(xy[1] - xy[0] - 1);
		ScalarFunc constFunc = (xy, label, param) -> 1.0;
		double[] intf = fs.integrate(new ScalarFunc[] {sourceFunc}, null);
		double[] mTh = fs.integrate(new ScalarFunc[] {constFunc}, null);
		System.out.println(intf[0] + "\t" + mTh[0]);
		ScalarFunc sourceFunc2 = (xy, label, param) -> sourceFunc.action(xy, label, param) - intf[0] / mTh[0];
		SRMatrix A = new SRMatrix(fs.GetNdof());
		double[] RHS = new double[fs.GetNdof()];
		fs.assembleStiff(new double[] {1.0}, BVPType.COMMON, A);
		fs.assembleSource(new ScalarFunc[] {sourceFunc2}, null, BVPType.COMMON, RHS);
		double[] x = new double[fs.GetNdof()];
		//A.sort();
		//System.out.println(A.display());
		fs.applyBC_Periodic(A, RHS, new double[] {1.0, -1.0}, 1, 3);
		fs.applyBC_Periodic(A, RHS, new double[] {1.0, 1.0}, 2, 4);
		//fs.applyBC_Periodic(A, RHS, new int[][] {{2, 4}, {1, 3}}, new Direct[] {Direct.X, Direct.Y});
//		System.out.println(A.display());
		//fs.applyBC_MBN(A, RHS, Direct.All, 0.0, 1, 3);
		//A.sort();
		//System.out.println(A.display());
		//System.out.println(A.isSymmetric());
		NewIterSSolver solver = new NewIterSSolver(A);
		solver.PCGSSOR(RHS, x, 1.5, 2);
//		DirectSSolver solver = new DirectSSolver(A, fs.getnBand());
//		NewIterSSolver solver = new NewIterSSolver(A);
//		solver.PLUSolve(RHS, x, 1);
		mesh2D.toTecplot("PeriodicSolution2.dat", x);		
	}

}
